Stability of Attractive Bose-Einstein Condensates in a Periodic Potential 
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Using a standing light wave trap, a stable quasi-one-dimensional attractive dilute-gas Bose- 
Einstein condensate can be realized. In a mean-field approximation, this phenomenon is modeled 
by the cubic nonlinear Schrodinger equation with attractive nonlinearity and an elliptic function 
potential of which a standing light wave is a special case. New families of stationary solutions are 
presented. Some of these solutions have neither an analog in the linear Schrodinger equation nor 
in the integrable nonlinear Schrodinger equation. Their stability is examined using analytic and 
numerical methods. Trivial-phase solutions are experimentally stable provided they have nodes and 
their density is localized in the troughs of the potential. Stable time-periodic solutions are also 
examined. 



I. INTRODUCTION 



(1) 



> 

&\ 
O 
(N 

o 
o 

-I— > 



o 
o 



X 



Dilute gas Bose-Einstein condensates (BECs) have 
been generated by many groups using different gases 
which arc cooled to very low temperatures and confined 
in magnetic fields or standing light waves. The sign of 
the atomic coupling determines whether the interaction 
of the BECs is repulsive or attractive. Note that ef- 
ficient tuning between attractive and repulsive conden- 
sates can be achieved via a Feshbach resonance Q. Re- 
pulsive BECs are experimentally stable j|. In contrast, 
attractive Lithium BECs have been shown to collapse in 
three dimensions |J-jq[, but are predicted to be stable in 
quasi-one-dimension |6j . Here, we study the dynamics 
and stability of quasi-one-dimcnsional, attractive BECs 
trapped in standing light waves. 

The mean-field description for the macroscopic BEC 
wavefunction is constructed using the Hartree-Fock ap- 
proximation [Q, resulting in the Gross-Pitaevskii equa- 
tion |^||. The quasi-one-dimensional regime of the 
Gross-Pitaevskii equation holds when the transverse di- 
mensions of the condensate are on the order of its heal- 
ing length and the longitudinal dimension is much longer 
than its transverse dimensions ||,[l0) . In this regime the 
BEC remains phase coherent and the governing equa- 
tions are one-dimensional. This is in contrast to a truly 
ID mean-field theory which requires transverse dimen- 
sions on the order of or less than the atomic interac- 
tion length pH| . In quasi-lD, the Gross-Pitaevskii equa- 
tion reduces to the cubic nonlinear Schrodinger equation 
(NLS) with a potential @,|||||. 

In this paper we construct new solutions correspond- 
ing to a quasi-one-dimensional attractive BEC trapped 
in an external periodic potential. The governing equa- 
tion is given by the nonlinear Schrodinger equation with 
a potential 



where ip(x,t) represents the macroscopic wave function 
of the condensate and V{x) is an experimentally gen- 
erated macroscopic potential. A large class of periodic 
potentials is given by 



V(x) = -V sn 2 (x,k) 



(2) 



where sn(x, k) denotes the Jacobian elliptic sine func- 
tion |l4| with elliptic modulus < k < 1. In the limit 
k — ► 1 — , V(x) becomes an array of well-separated hyper- 
bolic secant potential barriers or wells, while in the limit 
A: — ^ + it becomes purely sinusoidal. We note that for 
most intermediate values (e.g. k = 1/2) the potential 
closely resembles sinusoidal behavior and thus provides 
a good approximation to a standing wave potential. The 
potential is plotted in Fig. |l| for values of k = 0, 0.9, 0.999 
and 0.999999. Only for k very near unity (e.g. > 0.999) 
does the solution start to appear visibly elliptic. The free- 
dom in choosing k allows great flexibility in considering a 
wide variety of physically realizable periodic potentials. 

The paper is outlined as follows: in the next section we 
derive and consider various properties and limits of two 
types of explicit solutions of Eq. ([!]) with (0) . Section III 
develops the analytic framework for the linear stability 
properties of the new solutions of Sec. II. The stability 
results are confirmed by numerical computations. For 
many cases, the stability analysis yields only partial an- 
alytical results, and we rely on numerical experiments 
to determine stability. Nonstationary solutions are dis- 
cussed in Sec. IV and illustrated with various types of 
time-periodic solutions. We conclude the paper in Sec. V 
with a brief summary of the primary results of the paper 
and their consequences for the dynamics of an attractive 
BEC. 
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FIG. 1. The sn 2 (x,k) structure of the potential 
for varying values of k. Note that the ^-coordinate 
has been scaled by the period of the elliptic func- 
tion. Since sn(x, k) is periodic in x with period 

sin 2 a, V(x) is periodic in x 
with period 2K(k). This period approaches infinity as 
k — ► 1. 



II. STATIONARY SOLUTIONS 

Equation (Q) with V{x) — is an integrable equa- 
tion and many explicit solutions corresponding to vari- 
ous boundary conditions are known. A comprehensive 
overview of these solutions is found in If V{x) ^ 0, 
NLS is not integrable. In this case, only small classes of 
explicit solutions can most likely be obtained. Our choice 
of potential (||) is motivated by the form of the stationary 
solution of the NLS with V(x) — 0. An overview of these 
stationary solutions and their properties is found in JlO[ . 
At present, we restrict our attention to stationary solu- 
tions of Eq. (^), i.e., solutions whose time-dependence is 
restricted to 



tp(x,t) — r(x) exp(—iu>t -\- i6(x)) , 



(3) 



If 6 X = 0, then the solution is referred to as having trivial 
phase and we choose 9{x) = 0. Substituting the ansatz 
Eq. (H) in Eq. (|l|) and dividing out the exponential factor 
results in two equations: one from the real part and one 
from the imaginary part. The second equation can be 
integrated: 



(x) = c 



dx' 
o r 2 (x') ' 



(4) 



where c is a constant of integration. Note that 9(x) is a 
monotonous function of x. Substitution of this result in 
the remaining equation gives 



4 / \ c 2 r 3 (x)r"(x) 
^ {X) = J 2 



-r 6 (x) - V sn 2 (x,k)r 4 (x). 



(5) 



The following subsections describe two classes of solu- 
tions of this equation. 



Type A 

1. Derivation 

For these solutions, r 2 (x) is a quadratic function of 
sn(x, k): 



r 2 (x) = A sn 2 (x, k) + B. 



(6) 



Substituting this ansatz in Eq. (|5|) and equating the co- 
efficients of equal powers of sn(:r, k) results in relations 
among the solution parameters u>, c, A and B and the 
equation parameters Vq and k. These are 



W = \ I 1 



= B 



k 2 - 3B 
B 



BVn 



(7a) 



V + k 2 
A = -(V + k 2 ). 



V Q + k 2 

- lj (V Q + k 2 - Bk 2 ) , (7b) 

(7c) 



For a given potential V(x), this solution class has one 
free parameter B which plays the role of a constant back- 
ground level or offset. The freedom in choosing the po- 
tential gives a total of three free parameters: Vq, k and 
B. 

The requirements that both r 2 (x) and c 2 are positive 
imposes conditions on the domain of these parameters: 



Vq < -k 2 , B > 0, 
V > -k 2 , (V + k 2 ) < B < | I 



k 2 



(8a) 
(8b) 



The region of validity of these solutions is displayed in 
Fig. | 
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FIG. 2. The region of validity of the solutions of Type 
A is displayed shaded for a fixed value of k. The edges of 
these regions correspond to various trivial phase solutions. 
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FIG. 3. The amplitude of the trivial phase solutions of 
Type A versus the potential strength Vb- 



For typical values of Vo, fc and i3, the above equations 
give rise to solutions of Eq. (0) which are not periodic in 
x: r(x) is periodic with period 2K(k), whereas exp(i6(x)) 
is periodic with period T = 6^ 1 (2tt). In general these 
two periods 2K(k) and T are not commensurable. Thus, 
requiring periodic solutions results in another condition, 
namely 2K(k) JT — p/q, for two positive integers p and q. 
The most convenient way to express this phase quantiza- 
tion condition is to assume the potential (i.e., Vo and fc) is 
given, and to consider values of B for which the quantiza- 
tion condition is satisfied. Introducing j3 = BJ (Vb + fc 2 ), 
we find 



± 



vW-i)(i-fc 2 /3) r K w 



dx 



— sn(x, fc) 2 + f3 



P 

q 



(9) 



This equation is solved for (3, after which B = f3(Vo + k 2 ). 
For numerical simulations, the number of periods of the 
potential is set. This determines q, limiting the number 
of solutions of Eq. (^|) . Solutions with the same period- 
icity as the potential require p/q = 1. 

Note that solutions of Type A reduce to stationary so- 
lutions of Eqs. (P and (||) with V = 0. Furthermore, 
all stationary solutions of the integrable equation are ob- 
tained as limits of solutions of Type A. 



2. Limits and Properties 

The properties of these solutions are best understood 
by considering their various limit cases. 

The trivial phase case: The solutions of Type A 
have trivial phase when c = 0. Since c 2 has three fac- 
tors which are linear in B (see Eq. (R)), there are three 
choices of B for which this occurs: B — 0, B = Vb + fc 2 
and B = (Vb + fc 2 )/fc 2 . These possibilities are three of 
the four boundary lines of the region of validity in Fig. |^. 
Note that the remaining boundary line (Vb = — fc 2 ) cor- 
responds to r 2 (x) = B, which gives rise to a plane wave 
solution. Using Jacobian elliptic function identities [|l4j , 
one finds that the three other boundary lines give rise to 
simplified solution forms: B — gives 



ri(x) = \/-(V + k 2 ) sn(x, fc), oj 
B = Vq + fc 2 gives 



1 + fc 2 



(10) 



r-2 



(x) = VVo + k 2 cn(x,fc), w=i-V -fc 2 , (11) 



where cn(x, fc) denotes the Jacobian elliptic cosine func- 
tion. Lastly, B = (V + k 2 )/k 2 gives 



Wo + fc 2 
r 3\X) = : an{x,K), LO 



Vo 

fc 2 



fc 2 

2 ' 



(12) 



where dn(x, fc) denotes the third Jacobian elliptic func- 
tion. Solution ( |T~0| ) is valid for Vb < — fc 2 , whereas the 
other two solutions (11) and ( |l2] ) are valid for Vq > — fc 2 . 
The amplitude of these solutions as a function of poten- 
tial strength V is shown in Fig. ||. 

Both cn(x, fc) and sn(x, fc) have zero average as func- 
tions of x and lie in [—1, 1]. On the other hand, dn(a;, fc) 
has nonzero average. Its range is [VI — fc 2 , 1] . Further- 
more, cn(x, fc) and sn(ai, fc) are periodic in x with period 
AK{k), whereas dn(x, fc) is periodic with period 2K(k). 
Some solutions with trivial phase are shown in Fig. [|. 

The trigonometric limit: In the limit fc — * 0, the 
elliptic functions reduce to trigonometric functions and 
V{x) = -V sin 2 (x) = (V /2) cos(2a;) - V /2. Then 



r 2 (x) = -V sin 2 (a;) + B, w=--B. 
In this case, the phase integral Eq. (0) results in 



tan(0(x)) = ±Vi - v o/ B tan(x). 



(13) 



(14) 
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FIG. 4. Trivial phase solutions for k = 0.5. V(x) is 
indicated with a solid line. For the top figure Vo = — 1. For 
the bottom figure Vo = 1. 
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(a) V =-1,B=1 



(b) V =-1,B=0.1 




FIG. 5. Phase and amplitude of the trigonometric solu- 
tions. For all these figures, the solid line denotes V(x), the 
dashed line is r(x) and the dotted line is 8(x)/(2tt). Note 
that 6(x) becomes piecewise constant, as B approaches the 
boundary of the region of validity. Far away from this 
boundary, 0(x) is essentially linear. 



Note that this formula guarantees that the resulting so- 
lution is periodic with the same period as the potential, 
so no phase quantization is required. In the trigonomet- 
ric limit, the wedge between the two regions of validity 
in Fig. ^disappears. This is no surprise, as in this limit, 
dn(x, k) — > 1, and the third trivial phase solution reduces 
to a plane wave solution. The corner point of the region 
of validity also moves to the origin. Some trigonometric 
solutions are illustrated in Fig. ||. 

The solitary wave limit: k = 1. In this limit the 
elliptic functions reduce to hyperbolic functions. Specifi- 
cally, sn(x, k) — tanh(x). Hence in this limit, the poten- 
tial has only a single well or a single peak. Then Vo < 
gives rise to a repulsive potential, whereas Vo > gives 
rise to an attractive potential: V(x) = —Vq tanh 2 (x). 
In this case the phase 0(x) of Eq. ({|) can be calculated 
explicitly: 



r 2 (x) = -(V + 1) tanh 2 (x) + B, 
9(x) — \l — — — x+arctan 



V + l 



Vb + 1 



-B 



(15a) 

tanh(x) I , (15b) 



which is valid for B > and Vq < —1- This solution is a 
stationary solitary wave of depression on a positive back- 
ground. It is reminiscent of the gray soliton solution of 
the NLS equation with repulsive nonlinearity. Note that 
this solution exists with an attractive potential. One 
such solution is illustrated in Fig. ||a. Another solution 
exists when B = Vo + 1 > 0: r(x) = \/V~o + 1 sech(x) and 



8(x) — 0. This solution represents a stationary elevated 
solitary wave. It is a deformation of the bright soliton 
solution of the NLS equation with attractive nonlinear- 
ity. Depending on Vo, it exists in cither an attractive 
(— 1 < Vq < 0) or a repulsive (Vo > 0) potential. These 
solutions are shown in Fig. ||b-c. 

Understanding the solitary wave limit facilitates the 
understanding of what occurs for k — > 1. In this case 
the solutions of Type A reduce to a periodic train of soli- 
tons with exponentially small interactions as illustrated 
in Fig. ^|d. The exponentially small nature of the inter- 
actions can be exploited analytically as is done in Sec. IV 
for nonstationary solutions. 



Type B 



1. Derivation 



For these solutions, r 2 (x) is linear in sn(x, k) or 
dn(x, k). First we discuss the solution with sn(x, k). The 
quantities associated with this solution will be denoted 
with a subindex 1. The quantities associated with the 
dn(x, k) solution receive a subindex 2. 

Substituting 



r\(x) = ai sn(x, k) + bx , 



(16) 



in Eq. (|^) and equating different powers of sn(x, k) gives 
the relations: 



Vo 



UJl = 



(17a) 
(17b) 



^(16a 2 -fc 4 )(16a 2 -fc 2 ), (17c) 



_ 4a 2 



(17d) 



The class of potentials Eq. (f2|) is restricted by the first 
of these relations so that Vo is in the narrow range 
-3fc 2 /8 < Vo < 0. The solution class now depends on 
one free amplitude parameter ai and the free equation 
parameter k. 

The region of validity of this solution is, as before, de- 
termined by the requirements cf > and r\(x) > 0: 



(18) 
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(a) V =-2, B=0.1 



(b) V =-0.5, B=1 





(d) V =-2, B=1, k=1-10" 
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FIG. 6. Solutions with k=l (a,b,c) or k -» 1 (d). The 
solid line denotes V(x), the dashed line is r{x) and the dot- 
ted line is 6(x)/2tt. In (d), a value of k = 1 — 1CT 16 was 
used. 



The period of ri (x) is twice the period of the potential. 
Requiring periodicity in x of this first solution of Type B 
gives 



± 



v/(16a 2 -fc 4 )(fc 2 ~16a 2 ) [ 2K W dx 



Awk 3 



= P 



(19) 



For given fc and integers p, q, this equation is solved for 
ai. 

The dn(x, fc) solutions are found by substituting 



r%(x) = a 2 dn(x, fc) + b 2 , 



(20) 



in Eq. (||). Equating different powers of dn(a;, fc) imposes 
the following constraints on the parameters: 



(21a) 
(21b) 



V = 


-h 2 , 

8 ' 




UJ 2 = 


\ (1 + fc2) 


+ 6a|, 


4 = 


|(16a^- 


l)(16al + fc 2 


b 2 = 


-4a 2 .. 





(21c) 
(21d) 



The class of potentials (g) is restricted as for the previous 
solution by the first of these relations. The solution class 
again depends on one free amplitude parameter a 2 and 
the free equation parameter k. 

The region of validity of this solution is once more de- 
termined by the requirements c 2 > and r^ix) > 0: 



< a 2 < 



The period of r 2 (x) is equal to the period of the po- 
tential. Requiring periodicity in x of this second solution 
of Type B gives 



± 



v/(16a 2 -l)(16a 2 + fc 2 -l) [ K{k) dx 



4a 2 — dn(x, k) 



P 
(I 



(23) 



For given k and integers p, q, this equation is solved to 
determine a 2 . 

In contrast to solutions of Type A, solutions of Type 
B do not have a nontrivial trigonometric limit. In fact, 
for solutions of Type B, this limit is identical to the limit 
in which the potential strength Vb = — 3fc 2 /8 approaches 
zero. Thus it is clear that the solutions of Type B have 
no analogue in the integrable nonlinear Schrodinger equa- 
tion. However, other interesting limits do exist. 



2. Limits and Properties 



The trivial phase case: Trivial phase corresponds 
to c = 0. This occurs precisely at the boundaries of the 
regions of validity. For the first solution of Type B, there 
are four possibilities: ai — k 2 /4, a\ — k/4, a\ — — k 2 /A 
or a\ = —k/4. By replacing x by x+2K (k), one sees that 
the last two possibilities are completely equivalent to the 
first two, so only the first two need to be considered. For 
ai = k 2 /4, 



ri(x) 



-(1 + sn(x, k)) , lji 



fc 2 
4 ' 



(24) 



Equating a\ = k/4 gives 



-(1 + k sn(x. fc)) , lu\ 



1 1 



(25) 



For the second solution, there are four possibilities: a 2 = 
or a 2 — vl — fc 2 /4. The first one of these results in 
a zero solution. The second one gives an interesting 
trivial phase solutions. Let fc' = Vl k 2 . Then, for 
a 2 = VT^/4 = fc'/4, 



k' 



1 fc" 



(a) = — (dn(ar,fc) - k') , iv 2 = - + 



(26) 



(22) These solutions are shown in Fig. [?]. 
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FIG. 7. Solutions of Type B with trivial phase. The 
figures correspond to, from top to bottom, k — 0.5, k = 0.9 
and k = 0.999. The potential is indicated with a solid line. 
The other curves are: (1) |ri (a;) | with ai — k 2 /4, (2) ri(x) 
with a2 = k/i and (3) |t"2 (a;) | with = k' /4. 



0.5 - 









V(x) 


/ r(x) 

/ 

/ 

/ 






















5 


-2.5 2.5 

X 


t 



FIG. 8. Solitary wave solution of Type B. The poten- 
tial is indicated with a solid line. The dashed line solution 
corresponds to a — 0.3, the dotted line to a = —0.3. 



The hyperbolic limit: In this limit Vo = —3/8 and 
the potential is V(x) = 3 tanh(a;) 2 /8. The region of valid- 
ity for r 2 shrinks to a 2 = 0, giving only the zero solution. 
However, the region of validity for n shrinks to ax = 1/4, 
giving rise to a non-trivial shock- like solution with trivial 
phase: rf(x) — (1 + tanh(x))/4. This solution is shown 
in Fig. ||. As for the solutions of type A, the solitary 
wave limit gives an idea of the behavior of the solution 
for values of k — >• 1. 



III. STABILITY 

We have found a large number of new solutions to the 
governing Eqs. O) and (||). However, only solutions that 



are stable can be observed in experiments. In this sec- 
tion, we consider the stability of the different solutions. 
Both analytical and numerical results are presented for 
the solutions with trivial phase. In contrast, only numer- 
ical results are discussed for the nontrivial phase cases. 

The linear stability of the solution (|J) is investigated. 
To do so, the exact solutions are perturbed by: 

ip(x, t) = (r(x) + e<f>(x, t)) exp[i(0(a;) - tut)} (27) 

where e < 1 is a small parameter. Collecting terms at 
0(e) gives the linearized equation. Its real and imaginary 
parts are U = {U x , U 2 Y = {Re[<f>]Jm[<f)\f: 



U, 



,z.u = , [ ^ /_ 



u, 



(28) 



where 



L- 
S 



-d x - 



r(x) r{x) 



- 3r(x) 2 + V(x) - uj, (29a) 

- r{xf + V(x) - to, (29b) 

(29c) 



and J 



is a skew-symmetric matrix. The 



-1 

1 

operator L is Hermitian as are L + and L_ while S 
is anti-Hermitian. Considering solutions of the form 
U(.r, t) = U(x) exp(At) gives the eigenvalue problem 



CXJ = AU, 



(30) 



where C — JL and A is complex. If all A are imaginary, 
then linear stability is established. In contrast, if there 
is at least one eigenvalue with a positive real part, then 



instability results. Using the phase invariance tp 
of Eq. (Q), Noether's theorem |l6| gives 



C 





r(x) 



= 0, 



(31) 



which implies that L-r(x) = 0. Thus A = is in the 
spectrum of L_. For general solutions of the form (j^), 
determining the spectrum of the associated linearized 
eigenvalue problem ( p8|) is beyond the scope of current 
methods. However, some cases of trivial-phase solutions 
(c = 0) are amenable to analysis. 

The Hermitian operators L± are periodic Schrodinger 
operators and thus the spectra of these operators is real 
and consists of bands of continuous spectrum contained 
in [A±,oo) @. Here A± denote the ground state eigen- 
values of L± respectively. They are given by 



A- 



inf i 

1011=1 



(32) 



G 



where \\4>\\ 2 = (<f>\<f>)- From the relation L + = L-—2r(x) 2 
it follows that A + < A_. Also A_ < since A = is an 
eigenvalue of 

If A_ = 0, then L_ is non-negative and self-adjoint, 

i 

so we can define the non-negative square root, Li, via 

the spectral theorem [ |16| , and hence the Hermitian oper- 
i i — 

ator H = L^L+Lt_ can be constructed. The eigenvalue 
problem for C in Eq. ( ^0| ) is then equivalent to 

(H + \ 2 )<p = 0, (33) 

i 

with Ui — Lt<f- Denote the left-most point of the spec- 
trum of H by fio. If fio > then A 2 < and the eigenval- 
ues of C are imaginary and linear stability results. Since 

ill 
H = L^L+L^l and Li is positive, fj,Q > if and only if 

L + is non-negative. In contrast, if /io < then A 2 > 
and C has at least one pair of real eigenvalues with op- 
posite sign. This shows the existence of a growing mode 
leading to instability of the solution. 

This non-perturbative method distinguishes between 
two cases 

A. If r(x) > then Eq. (|^) implies r(x) is the ground 
state of L_ so that A_ = @ while A + < 0. Thus 
the solution (||) is unstable. 



B. 



If r(x) has a zero, it is no longer the ground 
state [|6[ and A_ < 0. Thus A_ and A + are both 
negative and the situation is indefinite. The non- 
perturbative methods are insufficient to determine 
linear stability or instability. 



While the solutions in Eqs. (12) and (|25f ) satisfy the in- 
stability criterion, the stability of the remaining solutions 
is undetermined by this method. We therefore rely on 
direct computations of the nonlinear governing Eqs. (Q) 
and (||) to determine stability. Where the numerics in- 
dicate stability, perturbative methods are used to con- 
firm this observation. For all computational simulations, 
twelve spatial periods are used. However, to better il- 
lustrate the dynamics, typically four spatial periods are 
plotted. Moreover, all computations are performed with 
white noise included in the initial data. 



A. trivial phase: Type A 
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FIG. 9. Unstable evolution of a Type A dn(x, k) solu- 
tion given by Eq. ( |l2|) over 40 time units with k — 0.7 and 
V = -0.3. 
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FIG. 10. Top view of the unstable evolution of Type A 
dn(a;, k) solutions given by Eq. (FF3) over 40 time units with 
k = 0.7 and with V = -0.3 (left), V = 0.0 (center), and 
Vo = 0.3 (right). 



instability occurs earlier. For the middle figure of Fig. [T(J, 
Vq = and we are solving the integrablc nonlinear 
Schrodinger equation. The instability of this solution 
is not surprising since dn(x, k) resembles a plane wave 
which is known to be modulationally unstable jOI. For 
increasing values of the amplitude A, the nonlinearity in 
(|lf) becomes more dominant, resulting in an earlier onset 
of the modulational instability. 



1. dn(x,k) 

In the case of the dn(x, k) solution (|T2|), r(x) > and 
the instability criterion A applies. The unstable evolu- 
tion is depicted in Figs. || and With Vq = —0.3 and 
k = 0.7 the onset of instability occurs at t ss 20. Fig- 
ure ^ depicts an overhead view of the dynamics with 
the onset of instability for k — 0.7 and for values of the 
potential strength Vq — —0.3,0,0.3. For increasing val- 
ues of the amplitude, A = y/Vo + k' 2 /k, the onset of 



2. sn(x,k) 

For the sn(x,k) solution ([To|), the indeterminate sta- 
bility criterion B results. Stability analysis for the linear 
Schrodinger equation suggests that this case is unstable 
since the density is localized on the peaks of the poten- 
tial (see Fig. ||). Figure |ll| illustrates this behavior and 
shows the onset of instability to occur for t w 20. 
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FIG. 11. Unstable evolution of a Type A sn(a:,fc) solu- 
tion given by Eq. ( [fo| ) over 40 time units with k — 0.7 and 
Vo = -1.0. 



3. cn(x, k) 

For the cn(x,k) solution ([TT|), the indeterminate sta- 
bility criterion B results once again. However, now the 
stability analysis for the linear Schrodinger equation sug- 
gests that two distinct cases must be considered. For 
Vq > 0, the density is localized on the peaks of the po- 
tential, and the solution is unstable as illustrated in the 
bottom of Fig. [l2[ The onset of instability occurs near 
t « 30. In contrast, for —k 2 < Vq < 0, the density is 
localized in the troughs of the potential which suggests 
that the solution might be stable. Indeed, as the top of 
Fig. |l2| illustrates, the cn(x, k) solution is stable in this 
regime. 




FIG. 12. Dynamics of Type A cn(:r, k) solutions given by 
Eq. (0) over 40 time units with k = 0.7 and with Vo = -0.3 
(top, stable) and Vo = 0.3 (bottom, unstable). 



B. Trivial Phase: Type B 

1. sn(x,k) 

For the Type B sn(a;, k) solution @) with a\ = fc 2 /4, 
the indeterminate stability criterion B applies. For this 
solution, the period of the density is twice the period of 
the potential so that the density is not localized in the 
troughs of the potential. Stability analysis for the linear 
Schrodinger equation suggests that this case is unstable 
as with the Type A sn(x, k) solution. Figure [l3] illus- 
trates this behavior and shows the onset of instability to 
occur for t » 200 for k = 0.5 and t » 15 for k = 0.999. 

The Type B sn(x,k) solution (|3) with a\ = fe/4 is 
nodeless. Thus r(x) > and the instability criterion A 
results. Figure [l4| illustrates this behavior and shows the 
onset of instability to occur at t w 20 for k = 0.5. 



2. dn{x,k) 

For the Type B dn(x,k) solution (|2(j|), the indetermi- 
nate stability criterion B results again. However, the sta- 
bility analysis for the linear Schrodinger equation again 
suggests that the solution might be stable since the den- 
sity is localized in the troughs of the potential. The dy- 
namics of this solution for k = 0.5 is illustrated in Fig. |l^. 

IV. NONSTATIONARY SOLUTIONS: SPATIALLY 
EXTENDED BREATHERS 

Equation (^) has many solutions describing conden- 
sates that oscillate in time. In this section, we construct 
such space and time periodic solutions in the large well 
separation limit (k — > 1). We consider only time-periodic 
solutions for which the condensate in each potential well 
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FIG. 13. Unstable evolution ol Type B sn(a;, k) solutions 
with ai = fe 2 /4 given by Eq. ( pi| ) over 300 time units for 
k = 0.5 (top) and over 20 time units for k = 0.999. 



oscillates with the same frequency. More time-periodic 
solutions will be considered elsewhere jlTj. 

Our solutions are obtained through a series of approx- 
imations. First, we assume the ansatz 



oo 

tp(x,t) = Aexp(-iu)t) sech(a: -&.(*)). (34) 

k— — oo 



This equation is motivated by the fact that a trivial-phase 
cn(x, k — > l)-solution can be written as y^-__ cp sech(x — 
AjK(k)) up to exponentially small terms describing the 
interaction between neighboring peaks. Note that the 
amplitude A and frequency u> in Eq. (34) are fixed across 
the condensate. 

Using the variational Lagrangian reduction ap- 
proach , effective equations for the motion of a single 
lump of condensate are derived. The dynamics of the 
center of a lump of condensate (t) are then given by 



FIG. 14. Unstable evolution of a Type B sn(:r, k) solu- 
tion with a± = fc/4 given by Eq. ( p5| ) over 40 time units 
with k — 0.5. 



0.05 
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FIG. 15. Stable evolution of the Type B dn(:r, k) solu- 
tion given by Eq. (|26[) over 40 time units with k — 0.5. 



the Newtonian equation of motion, = — dW(^k — 
£l)/d£,k where is the center position of the fcth po- 
tential well, with potential 



W(() = -V D( 2 (a-l3( 2 + 0(( 3 )) 



(35) 



where v is the average width of the condensate jT7| , a = 
4/15 and j3 — 2/63. For a small displacement the conden- 
sate undergoes near-harmonic oscillations with frequency 
y/-2aV v. 

After the self-interaction, the next largest contribu- 
tions are the nearest-neighbor interactions which arise 
from exponentially small tail overlaps. These interac- 
tions result in the following lattice differential equation 
(LDE) for the fcth position £fc(i), 
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FIG. 16. Stable vibrational mode of the condensate cor- 
responding to the fixed point P2 over 100 time units with 
Vo = —0.1. The initial conditions are a perturbation of the 
cnoidal wave solution ( |ll| ) with k = 0.9791048444 giving 
an initial lump-to-lump separation of = 6. 



-4A 3 (e 



1 ) + W'(£ k -e k ), (36) 



Here A£ k = £/c+i — 6c is the separation between centers of 
consecutive condensates. The same LDE can be derived 
from soliton perturbation theory ||l9| , ^o"f or a variational 
approach and corresponds to a Toda lattice |22] with 
additional on-site potentials due to W. 

We look for oscillatory solutions of Eq. (|3^) by consid- 
ering a Fourier expansion for £ k [p3[ : 



&(*) = ^2 a k{j)cos(jnt). 



(37) 



3=0 



We insert this ansatz into Eq. (|36j) and Taylor expand the 
exponentials, keeping terms to second order. Equating 
coefficients of cos(jTM), we obtain a recurrence relation 
between the amplitudes a k = a k (j = 1): 

a k+1 = 3/3 W 4 + (2 + 4a W - ti 2 ) a k - a fc _i , (38) 

where by using |A£°| = 2K(k) we find A = 
ftexp(Aft:(/c))/2A 2 and W ={V /32A 3 ) ey:p{2AK(k)). 

The recurrence relation Eq. ( |38| ) may be written as a 
two dimensional map Q (x k+ll y k+1 ) = F(x k ,y k ) by 
defining x k = a k and ?/ fe = au-i- 



x k - 
Vk- 



-!=3(3W x% 
-1 = £fc- 



(2 + 4a TFo - A 2 



y fc (39a) 
(39b) 



Fixed points of this map are calculated by solving 
(xo,2/o) = F{ x o,yo)- This results in three fixed points 
Pi = (xo,Vo) (t = 1,2,3): Pi = (0,0), P 2 = (x*,y*), and 
P3 = (—x*, — y*) where 



±1 



In 2 -4aW 
3/3 W : 



(40) 



FIG. 17. Stable two-period vibrational mode of the 
condensate corresponding to a period-two orbit over 100 
time units with Vo = —0.1. The initial conditions are 
a perturbation of the cnoidal wave solution ([n]) with 
k = 0.9791048444 giving an initial lump-to-lump separa- 
tion of A£fc = 6. 




FIG. 18. Stable evolution of a breathing mode of the 
condensate over 50 time units with Vo = —0.1. The initial 
conditions are a perturbation of the cnoidal wave solution 
(0) with k = 0.9791048444 giving an initial lump-to-lump 
separation of A£fc = 6. 



provided the root is real. From Eq. (37), fixed points P2 
and P3 correspond to vibrational modes in which each 
lump of condensate oscillates in time with the same am- 
plitude, frequency O, and phase. This dynamic is shown 
in Fig. Wn which is obtained from numerical simulation 
of Eq. (fin. 
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Period-two orbits of this map are calculated by 
solving (xo,yo) = F(F(xo,yoj) and are of the form 
{(x, y), (—x, —y)} 0]. This period two orbit corresponds 
to two alternating amplitudes (x and —x) of oscillation 
for consecutive condensates. Since the period two orbit 
is symmetric with respect to the origin, consecutive con- 
densates oscillate with the same magnitude but opposite 
phase. The dynamics for this case are show in Fig. |l7|. 

In a similar fashion, period-three and higher behavior 
can also be analyzed with the given ansatz Eq. (Q). In 
addition to allowing the lump positions to vary in time, 
we can also capture amplitude time variations. For sim- 
plicity, we illustrate the case where only the amplitudes 
vary in time. The construction of the appropriate LDE 
follows from previous methods. The fixed-point solution 
is illustrated in Fig. |l8| where lumps of condensate os- 
cillate in time with the same time-periodic amplitude, 
frequency, and phase. This solution type is referred to as 
an extended breather. Spatially localized breathers also 
exist and will be considered elsewhere fl7|] 



the potential, and adjacent density peaks are separated 
by nodes. Therefore, we have demonstrated within the 
mean field model the existence of at least one exper- 
imentally stable stationary state of an attractive BEC 
trapped in a standing light wave. In addition, some time- 
periodic condensates offer an alternative to experimental 
stabilization in a standing light wave potential. 
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V. SUMMARY AND CONCLUSIONS 

We considered the attractive nonlinear Schrodinger 
equation with an elliptic function potential as a model 
for a trapped, quasi-one-dimensional Bose-Einstein con- 
densate. Two new families of periodic solutions of this 
equation were found and their stability was investigated 
both analytically and numerically. Additionally, stable 
time-periodic solutions have been analyzed. 

Using perturbations with trivial phase (analysis) or 
perturbations with random phase (numerics), we find 
that stationary trivial-phase solutions are stable pro- 
vided they have nodes and their density is localized in 
the troughs of the potential. Nodeless solutions are un- 
stable with respect to this same class of perturbations. 
This is reminiscent of the modulational instability of the 
plane wave solution of the attractive integrablc nonlinear 
Schrodinger equation. 

Using random-phase perturbations, we find all 
nontrivial-phase solutions to be unstable. However, the 
time scale for the onset of instability for nontrivial-phase 
solutions varies significantly: nontrivial-phase solutions 
with parameter values close to values for trivial-phase 
solutions appear unaffected by the perturbation for long 
times. Other nontrivial-phase solutions go unstable more 
quickly. This implies that even stable trivial-phase solu- 
tions are structurally unstable: the smallest amount of 
phase ramping causes such solutions to lose their stabil- 
ity, albeit on time scale which may be longer than the 
lifetime of the BEC. 

This result implies self-focusing of any attractive sta- 
tionary condensate. However, the effects of self-focusing 
can be negligible on the lifetime of the BEC if there is 
no phase ramping, the density is localized in the wells of 
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